
module simulate

use grids
use params
use states
use searchnew
use utilityfs
use wage

implicit none

real(prec), allocatable :: inctypep(:),incp(:,:), conp(:,:), assp(:,:), savp(:,:), xvalp(:,:), hlth(:,:) !parent outcome records
real(prec), allocatable :: inck(:,:), conk(:,:), assk(:,:), savk(:,:), xvalk(:,:), labor(:,:) !kid outcome records
real(prec), allocatable :: xval_ins(:,:), xval_noins(:,:), xval_insPO(:,:) !values with and without insurance (for debugging welfare)
real(prec), allocatable :: ssip(:,:),ssik(:,:), COHp(:,:),COHk(:,:)
real(prec), allocatable :: alivemat(:,:),marstat(:,:),informal(:,:),muwgt(:,:),kappct(:,:) !other outcome records
real(prec), allocatable :: bequest(:,:),kterma(:,:) !terminal amounts
real(prec), allocatable :: visitwp(:,:), visitwk(:,:) !visits to different nodes
real(prec), allocatable :: maxass(:)
integer   , allocatable :: insd(:,:), insdPO(:,:)
real(prec), allocatable :: wstardiv(:,:),wstarmar(:,:)
integer	  , allocatable :: inctypek(:), simcohort(:)

real(prec), allocatable :: totvalp(:),totconp(:),totincp(:),totsavp(:),totSIp(:)
real(prec), allocatable :: totvalk(:),totconk(:),totinck(:),totsavk(:),totSIk(:)
real(prec), allocatable :: totmar(:),totmu(:), totkap(:)
real(prec), allocatable :: totins(:,:), totinsPO(:,:) !by parent income type
real(prec), allocatable :: totinf(:,:) !by kid income type
real(prec), allocatable :: totPT(:), totFT(:),totFT_sick(:),totFT_well(:),totPT_sick(:),totPT_well(:),totFTcohort(:,:)
real(prec), allocatable :: totbeq(:),totkTa(:)
real(prec), allocatable :: tothlthy(:), tothlthy_inck(:,:)
real(prec), allocatable :: totalive(:),totalivec(:,:),totalive_incp(:,:),totalive_inck(:,:)
integer, allocatable :: totcohort(:)

real(prec), allocatable :: avvalp(:),avconp(:),avincp(:),avsavp(:),avSIp(:)
real(prec), allocatable :: avvalk(:),avconk(:),avinck(:),avsavk(:),avSIk(:)
real(prec), allocatable :: avhlthy(:),avmar(:),avmu(:),avkap(:)
real(prec), allocatable :: avins(:),avinsPO(:) !by parent income type
real(prec), allocatable :: avinf(:)
real(prec) :: avinfoverall,avSIpoverall
real(prec) :: avFToverall,avFTwelloverall, avFTsickoverall
real(prec) :: avPToverall,avPTwelloverall, avPTsickoverall
real(prec) :: avwealthp50, avwealthk50
real(prec), allocatable :: avPT(:),avFT(:),avFT_sick(:),avFT_well(:),avPT_sick(:),avPT_well(:)
real(prec), allocatable :: wealthp10(:),wealthp25(:),wealthp50(:),wealthp75(:),wealthp90(:)
real(prec), allocatable :: wealthk10(:),wealthk25(:),wealthk50(:),wealthk75(:),wealthk90(:)
real(prec), allocatable :: wealthcp10(:,:),wealthcp25(:,:),wealthcp50(:,:),wealthcp75(:,:),wealthcp90(:,:)
real(prec), allocatable :: wealthck10(:,:),wealthck25(:,:),wealthck50(:,:),wealthck75(:,:),wealthck90(:,:)
real(prec), allocatable :: avFTcohort(:)
real(prec), allocatable :: kappap25(:),kappap50(:),kappap75(:)
real(prec), allocatable :: medconp(:),medconk(:)

real(prec), allocatable :: wlthmedian(:)

! integer, private :: seed

integer :: outcsv !output to csv files

contains

!------------------------------------------------------------------------------------------------------------------
subroutine allocate_simulate ( iret )
	implicit none
	integer, intent(out) :: iret
	integer		:: ier(1)
	allocate (maxass(nt),&
						inctypep(nosims),incp(nosims,nt),conp(nosims,nt),assp(nosims,nt),savp(nosims,nt),xvalp(nosims,nt),&
						xval_ins(nosims,nt),xval_noins(nosims,nt),xval_insPO(nosims,nt),&
						hlth(nosims,nt),insd(nosims,nt),insdPO(nosims,nt),wstardiv(nosims,nt),wstarmar(nosims,nt),&
						simcohort(nosims),&
						inctypek(nosims),inck(nosims,nt),conk(nosims,nt),assk(nosims,nt),savk(nosims,nt),xvalk(nosims,nt),labor(nosims,nt),&
						alivemat(nosims,nt),marstat(nosims,nt),informal(nosims,nt),muwgt(nosims,nt),&
						ssip(nosims,nt),ssik(nosims,nt),&
						COHp(nosims,nt),COHk(nosims,nt),&
						bequest(nosims,nt),kterma(nosims,nt), &
						kappct(nosims,nt),&
						visitwp(nwp,nt),visitwk(nw,nt),&
						totvalp(nt),totconp(nt),totincp(nt),totsavp(nt),totSIp(nt),&
						totvalk(nt),totconk(nt),totinck(nt),totsavk(nt),totSIk(nt),&
						totmar(nt),totmu(nt),totkap(nt),totbeq(nt),totkTa(nt),&
						tothlthy(nt),tothlthy_inck(nt,2),&
						totalive(nt),totalivec(nt,ncohort),totalive_incp(nt,nwp),totalive_inck(nt,2),&
						totcohort(ncohort),&
						totins(2,5),totinsPO(2,5),&
						totinf(2,5), &
						totPT(nt),totPT_sick(nt),totPT_well(nt),totFT(nt),totFT_sick(nt),totFT_well(nt),&
						totFTcohort(nt,ncohort),&
						avvalp(nt),avconp(nt),avincp(nt),avsavp(nt),avSIp(nt+1),&
						avvalk(nt),avconk(nt),avinck(nt),avsavk(nt),avSIk(nt+1),&
						avhlthy(nt),avmar(nt),avmu(nt),avkap(nt),&
						avPT(nt+1), avPT_well(nt+1),avPT_sick(nt+1),&
						avFT(nt+1), avFT_well(nt+1),avFT_sick(nt+1),&
						avins(6),avinsPO(6), avinf(6), &
						wealthp10(nt),wealthp25(nt),wealthp50(nt),wealthp75(nt),wealthp90(nt), &
						wealthk10(nt),wealthk25(nt),wealthk50(nt),wealthk75(nt),wealthk90(nt), &
						wealthcp10(ncohort,9),wealthcp25(ncohort,9),wealthcp50(ncohort,9), &
						wealthcp75(ncohort,9),wealthcp90(ncohort,9), &
						wealthck10(ncohort,9),wealthck25(ncohort,9),wealthck50(ncohort,9), &
						wealthck75(ncohort,9),wealthck90(ncohort,9), &						
						avFTcohort(ncohort),&
						kappap25(nt),kappap50(nt),kappap75(nt), &
						medconp(nt),medconk(nt),&
 			     	stat=ier(1))
		if ( ier(1).ne.0) then
			iret = -1
			return
		endif
		iret = 0
end subroutine allocate_simulate
!-------------------------------------------------------------------------------------------------------------------

!-------------------------------------------------------------------------------------------------------------------
subroutine simmodel

	implicit none

	include 'mpif.h'

	integer i,t,k !counter
	integer index
	
	!set variables used for simulation averages to zero
	totvalp = 0.d0
	totvalk = 0.d0
	totconp = 0.d0
	totconk = 0.d0
	totincp = 0.d0
	totinck = 0.d0
	totsavp = 0.d0
	totsavk = 0.d0
	totSIp  = 0.d0
	totSIk  = 0.d0
	tothlthy= 0.d0
	tothlthy_inck= 0.d0
	totinf  = 0.d0
	totins  = 0.d0
	totinsPO=0.d0
	totPT   = 0.d0
	totPT_sick = 0.d0
	totPT_well = 0.d0
	totFT   = 0.d0
	totFT_sick = 0.d0
	totFT_well = 0.d0
	totFTcohort = 0.d0
	totmar  = 0.d0
	totmu   = 0.d0
	totalive= 0.d0
	totalivec= 0.d0
	totalive_incp = 0.d0
	totalive_inck = 0.d0
	totcohort = 0
	totbeq  = 0.d0
	totkTa  = 0.d0
	totkap  = 0.d0

	avinf=0.d0
	avins=0.d0

	maxass=0d0
	xval_ins=0d0
	xval_noins=0d0
	xval_insPO=0d0
	inctypep=0d0
	incp=0d0
	conp=0d0
	assp=0d0
	savp=0d0
	xvalp=0d0
	hlth=0d0
	insd=0d0
	insdPO=0d0
	inctypek=0d0
	inck=0d0
	conk=0d0
	assk=0d0
	savk=0d0
	xvalk=0d0
	labor=0d0
	alivemat=0d0
	marstat=0d0
	informal=0d0
	muwgt=0d0
	ssip=0d0
	ssik=0d0
	COHp=0d0
	COHk=0d0
	bequest=0d0
	kterma=0d0
	kappct=0d0
	wstardiv=0d0
	wstarmar=0d0
	simcohort=0d0
	
	!***********************************
	!LOOP FOR DIFFERENT INDIVIDUALS
	!$OMP PARALLEL DO PRIVATE(index,mnode,mweight,t,k)
	do index=simnode1,simnode2 !1,nosims
		call run_simperson(index)
	enddo	!end loop individuals
  	!$OMP END PARALLEL DO
	! ***********************************
	
  	call MPI_ALLREDUCE(xval_ins,tempvecsim,recsizesim,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
  	xval_ins=reshape(tempvecsim,(/nosims,nt/))
  	call MPI_ALLREDUCE(xval_noins,tempvecsim,recsizesim,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
  	xval_noins=reshape(tempvecsim,(/nosims,nt/))
  	call MPI_ALLREDUCE(xval_insPO,tempvecsim,recsizesim,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
  	xval_insPO=reshape(tempvecsim,(/nosims,nt/))
  	call MPI_ALLREDUCE(xvalp,tempvecsim,recsizesim,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
  	xvalp=reshape(tempvecsim,(/nosims,nt/))
  	call MPI_ALLREDUCE(xvalk,tempvecsim,recsizesim,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
  	xvalk=reshape(tempvecsim,(/nosims,nt/))
  	call MPI_ALLREDUCE(conp,tempvecsim,recsizesim,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
  	conp=reshape(tempvecsim,(/nosims,nt/))
  	call MPI_ALLREDUCE(conk,tempvecsim,recsizesim,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
  	conk=reshape(tempvecsim,(/nosims,nt/))
  	call MPI_ALLREDUCE(incp,tempvecsim,recsizesim,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
  	incp=reshape(tempvecsim,(/nosims,nt/))
  	call MPI_ALLREDUCE(inck,tempvecsim,recsizesim,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
  	inck=reshape(tempvecsim,(/nosims,nt/))
  	call MPI_ALLREDUCE(savp,tempvecsim,recsizesim,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
  	savp=reshape(tempvecsim,(/nosims,nt/))
  	call MPI_ALLREDUCE(savk,tempvecsim,recsizesim,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
  	savk=reshape(tempvecsim,(/nosims,nt/))
  	call MPI_ALLREDUCE(ssip,tempvecsim,recsizesim,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
  	ssip=reshape(tempvecsim,(/nosims,nt/))
  	call MPI_ALLREDUCE(ssik,tempvecsim,recsizesim,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
  	ssik=reshape(tempvecsim,(/nosims,nt/))
  	call MPI_ALLREDUCE(labor,tempvecsim,recsizesim,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
  	labor=reshape(tempvecsim,(/nosims,nt/))
  	call MPI_ALLREDUCE(hlth,tempvecsim,recsizesim,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
  	hlth=reshape(tempvecsim,(/nosims,nt/))
  	call MPI_ALLREDUCE(marstat,tempvecsim,recsizesim,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
  	marstat=reshape(tempvecsim,(/nosims,nt/))
  	call MPI_ALLREDUCE(muwgt,tempvecsim,recsizesim,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
  	muwgt=reshape(tempvecsim,(/nosims,nt/))
  	call MPI_ALLREDUCE(bequest,tempvecsim,recsizesim,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
  	bequest=reshape(tempvecsim,(/nosims,nt/))
  	call MPI_ALLREDUCE(alivemat,tempvecsim,recsizesim,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
  	alivemat=reshape(tempvecsim,(/nosims,nt/))
  	call MPI_ALLREDUCE(insd,tempvecsimint,recsizesim,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,ierr)
  	insd=reshape(tempvecsimint,(/nosims,nt/))
  	call MPI_ALLREDUCE(insdPO,tempvecsimint,recsizesim,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,ierr)
  	insdPO=reshape(tempvecsimint,(/nosims,nt/))	
  	call MPI_ALLREDUCE(informal,tempvecsim,recsizesim,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
  	informal=reshape(tempvecsim,(/nosims,nt/))
  	call MPI_ALLREDUCE(assp,tempvecsim,recsizesim,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
  	assp=reshape(tempvecsim,(/nosims,nt/))
  	call MPI_ALLREDUCE(assk,tempvecsim,recsizesim,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
  	assk=reshape(tempvecsim,(/nosims,nt/))
  	call MPI_ALLREDUCE(kappct,tempvecsim,recsizesim,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
  	kappct=reshape(tempvecsim,(/nosims,nt/))
  	call MPI_ALLREDUCE(wstardiv,tempvecsim,recsizesim,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
  	wstardiv=reshape(tempvecsim,(/nosims,nt/))
  	call MPI_ALLREDUCE(wstarmar,tempvecsim,recsizesim,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
  	wstarmar=reshape(tempvecsim,(/nosims,nt/))
  	call MPI_ALLREDUCE(simcohort,tempvecsimint2,recsizesim2,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,ierr)
	simcohort=reshape(tempvecsimint2,(/nosims/))	
  	call MPI_ALLREDUCE(inctypek,tempvecsimint2,recsizesim2,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,ierr)
	inctypek=reshape(tempvecsimint2,(/nosims/))	

  do index=1,nosims
	do i=1,5
		if (simcohort(index)==i) then
		totcohort(i) = totcohort(i) + 1
		endif
	enddo
	do t=1,nt
  			!summing up for averages
			totvalp(t) = totvalp(t) + xvalp(index,t)
			totvalk(t) = totvalk(t) + xvalk(index,t)
			totconp(t) = totconp(t) + conp(index,t)
			totconk(t) = totconk(t) + conk(index,t)
			totincp(t) = totincp(t) + incp(index,t)
			totinck(t) = totinck(t) + inck(index,t)
			totsavp(t) = totsavp(t) + savp(index,t)
			totsavk(t) = totsavk(t) + savk(index,t)
			totSIp(t)  = totSIp(t)  + ssip(index,t)
			totSIk(t)  = totSIk(t)  + ssik(index,t)
			totkap(t)  = totkap(t)  + kappct(index,t)
			if (labor(index,t)==3 .and. alivemat(index,t)==1) then
				totFT(t) = totFT(t) + 1
				if (hlth(index,t).eq.1) then
					totFT_well(t)=totFT_well(t)+1
				else
					totFT_sick(t)=totFT_sick(t)+1
				endif
				do i=1,5
					if (simcohort(index)==i) then
						totFTcohort(t,i) = totFTcohort(t,i)+1
					endif
				enddo
			endif
			if (labor(index,t)==2 .and. alivemat(index,t)==1) then
				totPT(t) = totPT(t) + 1
				if (hlth(index,t).eq.1) then
					totPT_well(t)=totPT_well(t)+1
				else
					totPT_sick(t)=totPT_sick(t)+1
				endif
			endif
			if (hlth(index,t).eq.1) then
				tothlthy(t)= tothlthy(t) + 1
			endif
			totmar(t)  = totmar(t) + marstat(index,t)
			totmu(t)   = totmu(t) + muwgt(index,t)
			totbeq(t)  = totbeq(t) + bequest(index,t)
			totkTa(t)  = totkTa(t) + kterma(index,t)
			totalive(t)= totalive(t) + alivemat(index,t)
			do i=1,5
				if (simcohort(index)==i) then
					totalivec(t,i)= totalivec(t,i) + alivemat(index,t)
				endif
			enddo

			!informal care
			if (simcohort(index)==1) then
			if (assp(index,1)<1130.476 .and. (hlth(index,t).eq.2 .or. hlth(index,t).eq.3)) then
				totinf(1,1) = totinf(1,1) + informal(index,t)
				totinf(2,1) = totinf(2,1) + 1
			elseif (assp(index,1)<39336.48 .and. (hlth(index,t).eq.2 .or. hlth(index,t).eq.3)) then
				totinf(1,2) = totinf(1,2) + informal(index,t)
				totinf(2,2) = totinf(2,2) + 1
			elseif (assp(index,1)<118997.5 .and. (hlth(index,t).eq.2 .or. hlth(index,t).eq.3)) then
				totinf(1,3) = totinf(1,3) + informal(index,t)
				totinf(2,3) = totinf(2,3) + 1		
			elseif (assp(index,1)<296660.7 .and. (hlth(index,t).eq.2 .or. hlth(index,t).eq.3)) then
				totinf(1,4) = totinf(1,4) + informal(index,t)
				totinf(2,4) = totinf(2,4) + 1
			elseif (assp(index,1)>296660.7 .and. (hlth(index,t).eq.2 .or. hlth(index,t).eq.3)) then
				totinf(1,5) = totinf(1,5) + informal(index,t)
				totinf(2,5) = totinf(2,5) + 1
			endif	
			endif

		enddo

		!insurance by asset quintile
		if (simcohort(index)==1) then
			if (assp(index,1)<1130.476) then
				totins(1,1) = totins(1,1) + insd(index,1)
				totins(2,1) = totins(2,1) + 1
			elseif (assp(index,1)<39336.48) then
				totins(1,2) = totins(1,2) + insd(index,1)
				totins(2,2) = totins(2,2) + 1
			elseif (assp(index,1)<118997.5) then
				totins(1,3) = totins(1,3) + insd(index,1)
				totins(2,3) = totins(2,3) + 1			
			elseif (assp(index,1)<296660.7) then
				totins(1,4) = totins(1,4) + insd(index,1)
				totins(2,4) = totins(2,4) + 1
			elseif (assp(index,1)>296660.7) then
				totins(1,5) = totins(1,5) + insd(index,1)
				totins(2,5) = totins(2,5) + 1
			endif
		endif
		
		if (simcohort(index)==1) then
			if (assp(index,1)<1130.476) then
				totinsPO(1,1) = totinsPO(1,1) + insdPO(index,1)
				totinsPO(2,1) = totinsPO(2,1) + 1
			elseif (assp(index,1)<39336.48) then
				totinsPO(1,2) = totinsPO(1,2) + insdPO(index,1)
				totinsPO(2,2) = totinsPO(2,2) + 1
			elseif (assp(index,1)<118997.5) then
				totinsPO(1,3) = totinsPO(1,3) + insdPO(index,1)
				totinsPO(2,3) = totinsPO(2,3) + 1			
			elseif (assp(index,1)<296660.7) then
				totinsPO(1,4) = totinsPO(1,4) + insdPO(index,1)
				totinsPO(2,4) = totinsPO(2,4) + 1
			elseif (assp(index,1)>296660.7) then
				totinsPO(1,5) = totinsPO(1,5) + insdPO(index,1)
				totinsPO(2,5) = totinsPO(2,5) + 1
			endif
		endif

	enddo
	
	print*,'hi4'

if (DoEstimation==0 .and. pid==0) then
	!for every period, store max end-of-period savings in a vector for deciding upper bound of asset grid
	maxass = maxval(savp, DIM=1) !DIM=1 means columns (i.e.time) of the matrix sav
	open(unit=outcsv,file='output/simmaxassp.csv',action="write",status="replace")
	do i=1,nt
	    write(outcsv,"(100000(f0.12,',',:))") maxass(i)
	end do
	!for every period, store max end-of-period savings in a vector for deciding upper bound of asset grid
	maxass = maxval(savk, DIM=1) !DIM=1 means columns (i.e.time) of the matrix sav
	open(unit=outcsv,file='output/simmaxassk.csv',action="write",status="replace")
	do i=1,nt
	    write(outcsv,"(100000(f0.12,',',:))") maxass(i)
	end do

	open(unit=outcsv,file='output/simassp.csv',action="write",status="replace")
	do i=1,nosims
	    write(outcsv,"(100000(f0.12,',',:))") assp(i,:)
	enddo
	open(unit=outcsv,file='output/simassk.csv',action="write",status="replace")
	do i=1,nosims
	    write(outcsv,"(100000(f0.12,',',:))") assk(i,:)
	enddo
	open(unit=outcsv,file='output/simconp.csv',action="write",status="replace")
	do i=1,nosims
	    write(outcsv,"(100000(f0.12,',',:))") conp(i,:)
	enddo
	open(unit=outcsv,file='output/simconk.csv',action="write",status="replace")
	do i=1,nosims
	    write(outcsv,"(100000(f0.12,',',:))") conk(i,:)
	enddo
	open(unit=outcsv,file='output/simincp.csv',action="write",status="replace")
	do i=1,nosims
	    write(outcsv,"(100000(f0.12,',',:))") incp(i,:)
	enddo
	open(unit=outcsv,file='output/siminck.csv',action="write",status="replace")
	do i=1,nosims
	    write(outcsv,"(100000(f0.12,',',:))") inck(i,:)
	enddo
	open(unit=outcsv,file='output/simsavp.csv',action="write",status="replace")
	do i=1,nosims
	    write(outcsv,"(100000(f0.12,',',:))") savp(i,:)
	enddo
	open(unit=outcsv,file='output/simsavk.csv',action="write",status="replace")
	do i=1,nosims
	    write(outcsv,"(100000(f0.12,',',:))") savk(i,:)
	end do
	open(unit=outcsv,file='output/siminsd.csv',action="write",status="replace")
	do i=1,nosims
	    write(outcsv,"(100000(I0,',',:))") insd(i,:)
	end do
	open(unit=outcsv,file='output/siminsdPO.csv',action="write",status="replace")
	do i=1,nosims
	    write(outcsv,"(100000(I0,',',:))") insdPO(i,:)
	end do	
	open(unit=outcsv,file='output/simalive.csv',action="write",status="replace")
	do i=1,nosims
	    write(outcsv,"(100000(f0.12,',',:))") alivemat(i,:)
	end do
	open(unit=outcsv,file='output/simmarstat.csv',action="write",status="replace")
	do i=1,nosims
	    write(outcsv,"(100000(f0.12,',',:))") marstat(i,:)
	end do
	open(unit=outcsv,file='output/simmu.csv',action="write",status="replace")
	do i=1,nosims
	    write(outcsv,"(100000(f0.12,',',:))") muwgt(i,:)
	end do
	open(unit=outcsv,file='output/simkap.csv',action="write",status="replace")
	do i=1,nosims
	    write(outcsv,"(100000(f0.12,',',:))") kappct(i,:)
	end do
	open(unit=outcsv,file='output/simhlth.csv',action="write",status="replace")
	do i=1,nosims
	    write(outcsv,"(100000(f0.12,',',:))") hlth(i,:)
	end do
	open(unit=outcsv,file='output/simlabor.csv',action="write",status="replace")
	do i=1,nosims
	    write(outcsv,"(100000(f0.12,',',:))") labor(i,:)
	end do
	open(unit=outcsv,file='output/simfamily.csv',action="write",status="replace")
	do i=1,nosims
	    write(outcsv,"(100000(f0.12,',',:))") informal(i,:)
	end do
	open(unit=outcsv,file='output/simbeq.csv',action="write",status="replace")
	do i=1,nosims
	    write(outcsv,"(100000(f0.12,',',:))") bequest(i,:)
	end do
	open(unit=outcsv,file='output/simkterma.csv',action="write",status="replace")
	do i=1,nosims
	    write(outcsv,"(100000(f0.12,',',:))") kterma(i,:)
	end do

	open(unit=outcsv,file='output/simsocinsp.csv',action="write",status="replace")
	do i=1,nosims
	    write(outcsv,"(100000(f0.12,',',:))") ssip(i,:)
	end do

	open(unit=outcsv,file='output/simsocinsk.csv',action="write",status="replace")
	do i=1,nosims
	    write(outcsv,"(100000(f0.12,',',:))") ssik(i,:)
	end do

	open(unit=outcsv,file='output/simCOHp.csv',action="write",status="replace")
	do i=1,nosims
	    write(outcsv,"(100000(f0.12,',',:))") COHp(i,:)
	end do

	open(unit=outcsv,file='output/simCOHk.csv',action="write",status="replace")
	do i=1,nosims
	    write(outcsv,"(100000(f0.12,',',:))") COHk(i,:)
	end do

	open(unit=outcsv,file='output/simwstardiv.csv',action="write",status="replace")
	do i=1,nosims
	    write(outcsv,"(100000(f0.12,',',:))") wstardiv(i,:)
	end do
	open(unit=outcsv,file='output/simwstarmar.csv',action="write",status="replace")
	do i=1,nosims
	    write(outcsv,"(100000(f0.12,',',:))") wstarmar(i,:)
	end do
	
	open(unit=outcsv,file='output/simvalue_ins.csv',action="write",status="replace")
	do i=1,nosims
	    write(outcsv,"(100000(f0.12,',',:))") xval_ins(i,:)
	end do
	open(unit=outcsv,file='output/simvalue_noins.csv',action="write",status="replace")
	do i=1,nosims
	    write(outcsv,"(100000(f0.12,',',:))") xval_noins(i,:)
	end do	
	open(unit=outcsv,file='output/simvalue_insPO.csv',action="write",status="replace")
	do i=1,nosims
	    write(outcsv,"(100000(f0.12,',',:))") xval_insPO(i,:)
	end do	
	
	open(unit=outcsv,file='output/simcohort.csv',action="write",status="replace")
	do i=1,nosims
	    write(outcsv,"(100000(I0,',',:))") simcohort(i)
	end do
	
	open(unit=outcsv,file='output/siminctypek.csv',action="write",status="replace")
	do i=1,nosims
	    write(outcsv,"(100000(I0,',',:))") inctypek(i)
	end do
endif

end subroutine simmodel
!-------------------------------------------------------------------------------------------------------------------

! *******************************************************
subroutine run_simperson(index)

	implicit none

	integer HRSperson
	integer HRScohort
	integer cohortagestart
	integer incktype

	real(prec) COHpsim,COHksim
	real(prec) assetpsim,assetksim

	integer dead,diethispd,dielastpd
	integer index, node
	integer bnode,hnode

	real blah1,blah2,blah3,blah4


	integer :: k,t !counters

	!interpolation (on parent assets, kid assets, pareto weight)
	integer :: mina_p,mina_k,mina_w,minm
	integer :: nodep_round,nodek_round,nodem_round
	real(prec) :: lb_tri(3),ub_tri(3),xpts_tri(3)
	real(prec) :: lb_bi(2),ub_bi(2),xpts_bi(2)
	real(prec) :: wgt1,wgt2
	real(prec) :: apnode_interp,node_interpk,mnode_interp
	integer :: nodegridk
	integer lb_ap, lb_ak, ub_ap, ub_ak
	real(prec) :: x_ap, x_ak

	real(prec) :: xvaluep_divorced, xvaluek_divorced
	real(prec) :: xvaluep_married, xvaluek_married
	real(prec) :: LFPsimmat(ni), famsimmat(ni), marsimmat(ni), kapsimmat(ni)
	integer :: musimmat(ni)
	real(prec) :: xvalue,xvaluep,xvaluek !value
	real(prec) :: conexppsim,conexpksim !consumption
	real(prec) :: savsim,savpsim,savksim !savings at end of period
	real(prec) :: famsim,kapsim,LFPsim,marsim
	integer :: musim
	real(prec) :: beqsim,ktermasim
	integer :: musimstart
	integer    :: inssim
	integer    :: socinspsim,socinsksim
	real(prec) :: wstar_divorced, wstar_married

	integer :: wknode,wpnode

	real(prec) :: temp11,temp22,temp33

	real(prec), dimension(3,ni,na,nonodesk) :: VFpdiv, VFkdiv, Cpdiv, Ckdiv
	real(prec), dimension(6,na,nonodesk) :: VFpmar, VFkmar, VFmar, Cpmar, Ckmar, kapmar

	real(prec) :: v_adj,vp_adj,vk_adj,cp_adj,ck_adj,kap_adj,LFP_adj,fam_adj
	integer mu_adj

	real(prec) :: xkCVF_deadp(3)
	real(prec) :: xkCVF_divorced(3)
	real(prec) :: xCVF_married(6)
	real(prec) :: xvaluek_divorcedmat(ni), xvaluep_divorcedmat(ni)
	real(prec) :: xvaluekmat(ni), xvaluepmat(ni), xvaluemat(ni)
	real(prec) :: xconskmat(ni), xconspmat(ni)
	real(prec) :: LFPdiv(3)
	real(prec) :: LFPmar(6), fammar(6)
	integer :: locdivtemp(1), locdiv(1), locmar(1),locdivmat(ni)
	integer :: insdiv
	integer :: i1,i2, ins

	real(prec) :: xconsk, xconsp, xvalue_married
	integer simage, simageto


	
	!assign HRS person
	HRSperson=ceiling(dble(nHRS)*randprofile(index))		
	
	!start off alive
	dead=0
	diethispd=0
	!100% probability of death at end
	randsurv(index*(nt))=1.0d0

	!start off married
 	if (noncoop==1) then
 		marsim=0
 	else
		marsim=1
 	endif

	simcohort(index) = HRSprof_cohort(HRSperson)
	inctypek(index) = HRSprof_kinc(HRSperson)+1

	!***********************************
	!LOOP OVER LIFE CYCLE
	!-----------------------------------
	do t=1,(nt)

		simage = t
		simageto = t+1

		if (simage<HRSprof_age(HRSperson)) then
			goto 500
		endif
		
		!KEEPING TRACK OF NODES
		node = (index-1)*(nt) + t   ! To keep track of an individual's life (AGE*PERSON)

		diethispd=0

		!INCOME:
		!TRANSLATE RANDOM NUMBER FOR INCOME TO AN INCOME NODE
		if (simage.eq.1 .or. simage.eq.HRSprof_age(HRSperson)) then         
			!retrieve income for individual i at period 1 (parent)
			wpnode=HRSprof_pinc(HRSperson)	
			incktype=HRSprof_kinc(HRSperson)+1
			!find income draw for individual i at period 1 (kid)
			do k=(nw/2)*(incktype-1)+1,(nw/2)*incktype
				if (randinck(node).le.wuntot(k)) then
					wknode = k
					visitwk(wknode,simage) = visitwk(wknode,simage) + 1.d0
					goto 56
				endif
			enddo
			write (17,*) 'problem w/initial income: biggest wuntot, randinc, simage'
56		continue		
		else
			!find income draw for individual i at period t (kid)
			do k=(nw/2)*(incktype-1)+1,(nw/2)*incktype
				if (randinck(node).le.wcontotk(k,wknode,simage)) then
					wknode = k
					visitwk(wknode,simage) = visitwk(wknode,simage) + 1.d0
					goto 45
				endif
			enddo
			write (17,*) 'problem transit(w,wnode,simage): biggest wcontot, randinc, simage'
45		continue
		endif

		if (dead.eq.1) then
		!dead this period

			!died in a previous period
			alivemat(index,t) = 0

			!LOCATE ASSETS ON THE GRID
			assetksim = savksim*r				 
			!assets held (not including income yet!) are between mina and mina+1
			call find_ass(2, assetksim, simage, mina_k,wgt1,wgt2)
			!which income grid to use
			nodegridk = (wknode-1)*na
			!so: cash is equal to assets plus income, which is between(nodegrid+mina,nodegrid+mina+1)
			node_interpk=wgt1*mina_k+wgt2*(mina_k+1)

			!CALCULATE CONDITIONAL VALUE FUNCTIONS
			xkCVF_deadp(1)=interp_deadparent(simage,nodegridk,mina_k,node_interpk,valuekNT_deadparent) !1st=insurance,3rd=health (both placeholders, dont matter for kid)
			xkCVF_deadp(2)=interp_deadparent(simage,nodegridk,mina_k,node_interpk,valuekPT_deadparent) !1st=insurance,3rd=health (both placeholders, dont matter for kid)
			xkCVF_deadp(3)=interp_deadparent(simage,nodegridk,mina_k,node_interpk,valuekFT_deadparent) !1st=insurance,3rd=health (both placeholders, dont matter for kid)

			if (xkCVF_deadp(1)>=xkCVF_deadp(2) .and. xkCVF_deadp(1)>=xkCVF_deadp(3)) then
				xvaluek=xkCVF_deadp(1)
				conexpksim=interp_deadparent(simage,nodegridk,mina_k,node_interpk,conexpkNT_deadparent)
				LFPsim=1
			else if (xkCVF_deadp(2)>=xkCVF_deadp(1) .and. xkCVF_deadp(2)>=xkCVF_deadp(3)) then
				xvaluek=xkCVF_deadp(2)
				conexpksim=interp_deadparent(simage,nodegridk,mina_k,node_interpk,conexpkPT_deadparent)
				LFPsim=2
			else
				xvaluek=xkCVF_deadp(3)
				conexpksim=interp_deadparent(simage,nodegridk,mina_k,node_interpk,conexpkFT_deadparent)
				LFPsim=3
			endif

			!SAVINGS FOR NEXT PERIOD
			if (LFPsim.eq.3) then
				COHksim = max(assetksim+incomek(wknode,simage) +nonwageinc(wknode)          ,cmink)
			elseif (LFPsim.eq.2) then
				COHksim = max(assetksim+PTwagefrac*incomek(wknode,simage)+nonwageinc(wknode),cmink)
			else
				COHksim = max(assetksim                                  +nonwageinc(wknode),cmink)		
			endif			

			if (COHksim==cmink) then
				socinsksim=1
			else
				socinsksim=0
			endif

			savksim=COHksim-conexpksim
			if (savksim<0) then
				conexpksim=conexpksim+savksim
				savksim=0
			endif

			!record matrix of all sorts of values
			!inctypek(index)  = incktype
			xvalk(index,t)   = xvaluek
			conk(index,t)    = conexpksim
			savk(index,t)	 = savksim
			assk(index,t)	 = assetksim
			COHk(index,t)	 = COHksim
			ssik(index,t)    = socinsksim
			inck(index,t)    = incomek(wknode,simage)
			labor(index,t)   = LFPsim
			
		else
		!alive this period

			!HEALTH: translate random number for health to health state
			if (simage-HRSprof_age(HRSperson)+1<=ntHRS) then
				hnode=HRSprof_health(HRSperson,simage-HRSprof_age(HRSperson)+1)
			else
				if (randhlth(node).le.hprob(1,hnode,wpnode,simage)) then
					hnode=1 !healthy this period
				else if (randhlth(node).gt.hprob(1,hnode,wpnode,simage) .and. &
							   randhlth(node).le.(hprob(1,hnode,wpnode,simage)+hprob(2,hnode,wpnode,simage))) then
					hnode=2 !semi-unhealthy this period
				else
					hnode=3 !unhealthy this period
				endif
			endif	

			!so big node:
			bnode = (hnode-1)*nwp + wpnode

			!not dead yet
			alivemat(index,t) = 1

			!check if person dies at end of this period
			if (simage == nt) then
				diethispd = 1 
			else if (simage-HRSprof_age(HRSperson)+1 <=ntHRS-1) then
			 if (HRSprof_health(HRSperson,simage-HRSprof_age(HRSperson)+1+1) == 0) diethispd = 1
			else if (simage-HRSprof_age(HRSperson)+1 >=ntHRS .and. randsurv(node)>=survprob(hnode,wpnode,simage)) then
				diethispd = 1
			endif

			!FIGURE OUT ASSETS AND PARETO WEIGHT: 
			!find initial assets, and set assets at start of each subsequent period
			!find initial pareto weight, and set pareto weight at start of each subsequent period
			if (simage.eq.1 .or. simage.eq.HRSprof_age(HRSperson)) then

				!parent assigned HRS person profile
				assetpsim=HRSprof_passets(HRSperson)
				if (assetpsim==0) assetpsim=1
				!kid also assigned HRS person profile (imputed assets)
				assetksim=HRSprof_kassets(HRSperson)
				if (assetksim==0) assetksim=1
				
				musimstart=musimstartparam
				! for sensitivity of weight to starting at different ages:
				!if (simage.eq.1) then
				!	musimstart=musimstartparam
				!else
				!	musimstart=musimstartparamc(simcohort(index)-1)
				!endif

			else    ! simage.gt.1
				assetpsim = savpsim*r
				assetksim = savksim*r				
				musimstart = musim !start with mu from previous period
			endif

			!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
			!LOCATE CONTINUOUS STATE VARIABLES ON THE GRID

			!parent assets:
			!assets held are between mina and mina+1
			call find_ass(1, assetpsim, simage, mina_p,wgt1,wgt2)
	    	apnode_interp=wgt1*mina_p+wgt2*(mina_p+1)
	    	lb_ap=mina_p
	    	ub_ap=mina_p+1
	    	x_ap=apnode_interp

			!kid assets:
			!assets held (not including income yet!) are between mina and mina+1
			call find_ass(2, assetksim, simage, mina_k,wgt1,wgt2)
			!which income grid to use
			nodegridk = (wknode-1)*na
			!so: cash is equal to assets plus income, which is between(nodegrid+mina,nodegrid+mina+1)
	    	node_interpk=wgt1*mina_k+wgt2*(mina_k+1)
	    	lb_ak=nodegridk+mina_k
	    	ub_ak=nodegridk+mina_k+1
	    	x_ak=nodegridk+node_interpk
			!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

			!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
			!CALCULATE CONDITIONAL VALUE FUNCTIONS
			!this depends on insurance, which is a state variable for simage>1 but choice for simage==1

			if (simage==1) then
				i1=1
				i2=ni
			elseif (simage>1 .and. simage.eq.HRSprof_age(HRSperson)) then
				!these are later cohorts so ltci is pre-determined
				i1=2-HRSprof_ltci(HRSperson)
				i2=2-HRSprof_ltci(HRSperson)
			else
				i1=inssim
				i2=inssim
			endif			

			lb_bi=(/lb_ap,lb_ak/)
			ub_bi=(/ub_ap,ub_ak/)
			xpts_bi=(/x_ap,x_ak/)

			!DIVORCE: 
			xvaluek_divorcedmat=-9999999
			xvaluep_divorcedmat=-9999999

			do ins=i1,i2

				!1. calculate conditional values for kid
				!   involves bilinear interpolation on parent assets, kid assets

				LFPdiv=(/1,2,3/)

				VFpdiv(1,ins,:,:)=valuepNT_divorced(simage,ins,bnode,:,:)
				VFpdiv(2,ins,:,:)=valuepPT_divorced(simage,ins,bnode,:,:)
				VFpdiv(3,ins,:,:)=valuepFT_divorced(simage,ins,bnode,:,:)

				VFkdiv(1,ins,:,:)=valuekNT_divorced(simage,ins,bnode,:,:)
				VFkdiv(2,ins,:,:)=valuekPT_divorced(simage,ins,bnode,:,:)
				VFkdiv(3,ins,:,:)=valuekFT_divorced(simage,ins,bnode,:,:)

				Cpdiv(1,ins,:,:)=conexppNT_divorced(simage,ins,bnode,:,:)
				Cpdiv(2,ins,:,:)=conexppPT_divorced(simage,ins,bnode,:,:)
				Cpdiv(3,ins,:,:)=conexppFT_divorced(simage,ins,bnode,:,:)

				Ckdiv(1,ins,:,:)=conexpkNT_divorced(simage,ins,bnode,:,:)
				Ckdiv(2,ins,:,:)=conexpkPT_divorced(simage,ins,bnode,:,:)
				Ckdiv(3,ins,:,:)=conexpkFT_divorced(simage,ins,bnode,:,:)

				xkCVF_divorced(1)=bilinear(xpts_bi,lb_bi,ub_bi, &
												  VFkdiv(1,ins,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
				xkCVF_divorced(2)=bilinear(xpts_bi,lb_bi,ub_bi, &
													VFkdiv(2,ins,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
				xkCVF_divorced(3)=bilinear(xpts_bi,lb_bi,ub_bi, &
													VFkdiv(3,ins,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))

				!2. kid chooses which LFP, use that conditional value function for both parent and child
				xvaluek_divorcedmat(ins)=maxval(xkCVF_divorced)
				locdivtemp=maxloc(xkCVF_divorced)
				locdivmat(ins)=locdivtemp(1)
				
				!3. parent calculates value based on kid's LFP
				xvaluep_divorcedmat(ins)=bilinear(xpts_bi,lb_bi,ub_bi, &
																 VFpdiv(locdivmat(ins),ins,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))

			enddo

			if (simage==1 .and. hnode.ne.1) then
					xvaluep_divorced=xvaluep_divorcedmat(2)
					xvaluek_divorced=xvaluek_divorcedmat(2)
					locdiv=locdivmat(2)
					insdiv=2	
			else  
				if (ni==3) then
					if (xvaluep_divorcedmat(1)>=xvaluep_divorcedmat(2) .and. xvaluep_divorcedmat(1)>=xvaluep_divorcedmat(3)) then
						xvaluep_divorced=xvaluep_divorcedmat(1)
						xvaluek_divorced=xvaluek_divorcedmat(1)
						locdiv=locdivmat(1)
						insdiv=1 !this should not change from previous period if simage>1...
					elseif (xvaluep_divorcedmat(3)>=xvaluep_divorcedmat(1) .and. xvaluep_divorcedmat(3)>=xvaluep_divorcedmat(2)) then
						xvaluep_divorced=xvaluep_divorcedmat(3)
						xvaluek_divorced=xvaluek_divorcedmat(3)
						locdiv=locdivmat(3)
						insdiv=3 !this should not change from previous period if simage>1...
					else
						xvaluep_divorced=xvaluep_divorcedmat(2)
						xvaluek_divorced=xvaluek_divorcedmat(2)
						locdiv=locdivmat(2)
						insdiv=2 !this should not change from previous period if simage>1...
					endif
				else
					if (xvaluep_divorcedmat(1)>=xvaluep_divorcedmat(2)) then
						xvaluep_divorced=xvaluep_divorcedmat(1)
						xvaluek_divorced=xvaluek_divorcedmat(1)
						locdiv=locdivmat(1)
						insdiv=1 !this should not change from previous period if simage>1...
					else
						xvaluep_divorced=xvaluep_divorcedmat(2)
						xvaluek_divorced=xvaluek_divorcedmat(2)
						locdiv=locdivmat(2)
						insdiv=2 !this should not change from previous period if simage>1...
					endif
				endif

				!HERE DO WTP CALCULATIONS FOR NON-COOPERATION COUNTERFACTUAL
				!find wealth value w* s.t. xvaluep_divorcedmat(1)=xvaluep_divorcedmat(2,w*)
				call find_wstar_div(xvaluep_divorcedmat(1),mina_w,wgt1,wgt2)
				wstar_divorced=wgt1*mina_w+wgt2*(mina_w+1)

			endif
			!done with divorce

			if (marsim==0) then
				xvaluep=xvaluep_divorced
				xvaluek=xvaluek_divorced
				conexppsim=bilinear(xpts_bi,lb_bi,ub_bi,Cpdiv(locdiv(1),insdiv,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
				conexpksim=bilinear(xpts_bi,lb_bi,ub_bi,Ckdiv(locdiv(1),insdiv,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
				inssim=insdiv
				famsim=0
				marsim=0
			
			else

			!MARRIAGE:
			xvaluekmat=-9999999
			xvaluepmat=-9999999
			xvaluemat=-9999999
			marsimmat=-9999
			do ins=i1,i2
				!1. calculate conditional values for kid and parent and overall
				!   involves trilinear interpolation on parent assets, kid assets, pareto weight

				LFPmar=(/1,1,2,2,3,3/)
				fammar=(/0,1,0,1,0,1/)

				VFpmar(1,:,:)=valuepNT0_married(simage,ins,bnode,musimstart,:,:)
				VFpmar(2,:,:)=valuepNT100_married(simage,ins,bnode,musimstart,:,:)
				VFpmar(3,:,:)=valuepPT0_married(simage,ins,bnode,musimstart,:,:)
				VFpmar(4,:,:)=valuepPT100_married(simage,ins,bnode,musimstart,:,:)
				VFpmar(5,:,:)=valuepFT0_married(simage,ins,bnode,musimstart,:,:)
				VFpmar(6,:,:)=valuepFT100_married(simage,ins,bnode,musimstart,:,:)

				VFkmar(1,:,:)=valuekNT0_married(simage,ins,bnode,musimstart,:,:)
				VFkmar(2,:,:)=valuekNT100_married(simage,ins,bnode,musimstart,:,:)
				VFkmar(3,:,:)=valuekPT0_married(simage,ins,bnode,musimstart,:,:)
				VFkmar(4,:,:)=valuekPT100_married(simage,ins,bnode,musimstart,:,:)
				VFkmar(5,:,:)=valuekFT0_married(simage,ins,bnode,musimstart,:,:)
				VFkmar(6,:,:)=valuekFT100_married(simage,ins,bnode,musimstart,:,:)

				VFmar(1,:,:)=valueNT0_married(simage,ins,bnode,musimstart,:,:)
				VFmar(2,:,:)=valueNT100_married(simage,ins,bnode,musimstart,:,:)
				VFmar(3,:,:)=valuePT0_married(simage,ins,bnode,musimstart,:,:)
				VFmar(4,:,:)=valuePT100_married(simage,ins,bnode,musimstart,:,:)
				VFmar(5,:,:)=valueFT0_married(simage,ins,bnode,musimstart,:,:)
				VFmar(6,:,:)=valueFT100_married(simage,ins,bnode,musimstart,:,:)

				Cpmar(1,:,:)=conexppNT0_married(simage,ins,bnode,musimstart,:,:)
				Cpmar(2,:,:)=conexppNT100_married(simage,ins,bnode,musimstart,:,:)
				Cpmar(3,:,:)=conexppPT0_married(simage,ins,bnode,musimstart,:,:)
				Cpmar(4,:,:)=conexppPT100_married(simage,ins,bnode,musimstart,:,:)
				Cpmar(5,:,:)=conexppFT0_married(simage,ins,bnode,musimstart,:,:)
				Cpmar(6,:,:)=conexppFT100_married(simage,ins,bnode,musimstart,:,:)
				
				Ckmar(1,:,:)=conexpkNT0_married(simage,ins,bnode,musimstart,:,:)
				Ckmar(2,:,:)=conexpkNT100_married(simage,ins,bnode,musimstart,:,:)
				Ckmar(3,:,:)=conexpkPT0_married(simage,ins,bnode,musimstart,:,:)
				Ckmar(4,:,:)=conexpkPT100_married(simage,ins,bnode,musimstart,:,:)
				Ckmar(5,:,:)=conexpkFT0_married(simage,ins,bnode,musimstart,:,:)
				Ckmar(6,:,:)=conexpkFT100_married(simage,ins,bnode,musimstart,:,:)

				kapmar(1,:,:)=kappaNT0_married(simage,ins,bnode,musimstart,:,:)
				kapmar(2,:,:)=kappaNT100_married(simage,ins,bnode,musimstart,:,:)
				kapmar(3,:,:)=kappaPT0_married(simage,ins,bnode,musimstart,:,:)
				kapmar(4,:,:)=kappaPT100_married(simage,ins,bnode,musimstart,:,:)
				kapmar(5,:,:)=kappaFT0_married(simage,ins,bnode,musimstart,:,:)
				kapmar(6,:,:)=kappaFT100_married(simage,ins,bnode,musimstart,:,:)

				xCVF_married(1)=bilinear(xpts_bi,lb_bi,ub_bi,VFmar(1,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
				xCVF_married(2)=bilinear(xpts_bi,lb_bi,ub_bi,VFmar(2,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
				xCVF_married(3)=bilinear(xpts_bi,lb_bi,ub_bi,VFmar(3,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
				xCVF_married(4)=bilinear(xpts_bi,lb_bi,ub_bi,VFmar(4,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
				xCVF_married(5)=bilinear(xpts_bi,lb_bi,ub_bi,VFmar(5,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
				xCVF_married(6)=bilinear(xpts_bi,lb_bi,ub_bi,VFmar(6,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))


				!2. choose max conditional value function
				xvalue_married=maxval(xCVF_married)
				locmar=maxloc(xCVF_married)
				xvaluep_married=bilinear(xpts_bi,lb_bi,ub_bi,VFpmar(locmar(1),int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
				xvaluek_married=bilinear(xpts_bi,lb_bi,ub_bi,VFkmar(locmar(1),int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))

				!CHECK PARTICIPATION CONSTRAINTS (4 SCENARIOS):

				!1. both prefer marriage:
				if (limcom==0 .or. (xvaluep_married >= xvaluep_divorced .and. xvaluek_married >= xvaluek_divorced)) then
					xvaluemat(ins)=xvalue_married
					xvaluepmat(ins)=xvaluep_married
					xvaluekmat(ins)=xvaluek_married
					xconspmat(ins)=bilinear(xpts_bi,lb_bi,ub_bi,Cpmar(locmar(1),int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
					xconskmat(ins)=bilinear(xpts_bi,lb_bi,ub_bi,Ckmar(locmar(1),int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
					LFPsimmat(ins)=LFPmar(locmar(1))
					famsimmat(ins)=fammar(locmar(1))
					marsimmat(ins)=1
					kapsimmat(ins)=bilinear(xpts_bi,lb_bi,ub_bi,kapmar(locmar(1),int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
					musimmat(ins)=musimstart !not modified

				!2.kid prefers divorce
				else if (xvaluek_married < xvaluek_divorced) then
					call find_valkwgt_sim(xvaluek_divorced, &
										  v_adj,vp_adj,vk_adj,cp_adj,ck_adj,kap_adj,mu_adj,LFP_adj,fam_adj)
					!remain married with modified pareto weight
					xvaluemat(ins)=v_adj
					xvaluepmat(ins)=vp_adj
					xvaluekmat(ins)=vk_adj
					xconspmat(ins)=cp_adj
					xconskmat(ins)=ck_adj
					LFPsimmat(ins)=LFP_adj
					famsimmat(ins)=fam_adj
					musimmat(ins)=mu_adj
					kapsimmat(ins)=kap_adj
					marsimmat(ins)=1
				
				!3.parent prefers divorce
				else
					call find_valpwgt_sim(xvaluep_divorced, &
										  v_adj,vp_adj,vk_adj,cp_adj,ck_adj,kap_adj,mu_adj,LFP_adj,fam_adj)
					!remain married with modified pareto weight
					xvaluemat(ins)=v_adj
					xvaluepmat(ins)=vp_adj
					xvaluekmat(ins)=vk_adj
					xconspmat(ins)=cp_adj
					xconskmat(ins)=ck_adj
					LFPsimmat(ins)=LFP_adj
					famsimmat(ins)=fam_adj
					musimmat(ins)=mu_adj
					kapsimmat(ins)=kap_adj
					marsimmat(ins)=1
				endif

			enddo !over insurance (choices if simage==1, state if simage>1)

			if (simage==1 .and. hnode.ne.1) then
				inssim=2
				xvalue=xvaluemat(inssim)
				xvaluep=xvaluepmat(inssim)
				xvaluek=xvaluekmat(inssim)
				conexppsim=xconspmat(inssim)
				conexpksim=xconskmat(inssim)
				LFPsim=LFPsimmat(inssim)
				famsim=famsimmat(inssim)
				if (marsimmat(2)==1) then
					musim=musimmat(inssim)
					kapsim=kapsimmat(inssim)
					marsim=1  		
				else 
					famsim=0
					marsim=0
				endif
			else
				if (marsimmat(1)==1) then
					temp11=paretoweight(musimmat(1))*xvaluepmat(1) + &
						  (1.0-paretoweight(musimmat(1)))*xvaluekmat(1)
				else
					temp11=-999999
			    endif
				if (marsimmat(2)==1) then
					temp22=paretoweight(musimmat(2))*xvaluepmat(2) + &
						  (1.0-paretoweight(musimmat(2)))*xvaluekmat(2)	
				else
					temp22=-999999
				endif
				if (marsimmat(1)==1 .and. marsimmat(2).ne.1) then
					inssim=1
				elseif (marsimmat(2)==1 .and. marsimmat(1).ne.1) then
					inssim=2
				elseif (marsimmat(1)==1 .and. marsimmat(2)==1) then
					if (temp11>=temp22) then
						inssim=1
					else
						inssim=2
					endif
				endif
				temp33=-999999
				if (ni==3) then
					if (marsimmat(3)==1) then
						temp33=paretoweight(musimmat(3))*xvaluepmat(3) + &
							(1.0-paretoweight(musimmat(3)))*xvaluekmat(3)
						if (inssim==1) then
							if (temp33>=temp11) then
								inssim=3
							endif
						elseif (inssim==2) then
							if (temp33>=temp22) then
								inssim=3
							endif
						else
							inssim=3
						endif
					endif
			    endif
				xvalue=xvaluemat(inssim)
				xvaluep=xvaluepmat(inssim)
				xvaluek=xvaluekmat(inssim)
				conexppsim=xconspmat(inssim)
				conexpksim=xconskmat(inssim)
				LFPsim=LFPsimmat(inssim)
				famsim=famsimmat(inssim)
				kapsim=kapsimmat(inssim)
				musim=musimmat(inssim)
				marsim=1
				if (marsimmat(1).ne.1 .and. marsimmat(2).ne.1) then
					if (ni==3) then
						if (marsimmat(3).ne.1) then
							inssim=insdiv
							xvalue=xvaluemat(insdiv)
							xvaluep=xvaluepmat(insdiv)
							xvaluek=xvaluekmat(insdiv)
							conexppsim=xconspmat(insdiv)
							conexpksim=xconskmat(insdiv)
							LFPsim=LFPdiv(locdiv(1))
							famsim=0
							marsim=0						
						endif
					else
						print*,"here 3"
						inssim=insdiv
						xvalue=xvaluemat(insdiv)
						xvaluep=xvaluepmat(insdiv)
						xvaluek=xvaluekmat(insdiv)
						conexppsim=xconspmat(insdiv)
						conexpksim=xconskmat(insdiv)
						LFPsim=LFPdiv(locdiv(1))
						famsim=0
						marsim=0
					endif
				endif
			endif

				if (simage==1) then
					!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
					!WTP CALCULATIONS FOR NON-COOPERATION COUNTERFACTUAL
					!find wealth value w* s.t. xvaluemat(1)=xvaluemat(2,w*) (not just for parent, but overall value!)

					!give money to parent:
					call find_wstar_mar(1,xvaluemat(1),mina_w,wgt1,wgt2)
					wstar_married=wgt1*mina_w+wgt2*(mina_w+1)
				
					!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
				endif

			endif !whether skip marriage stuff or not
			!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

			if (inssim.eq.1) then
				COHpsim = max(assetpsim+incomep(wpnode,simage)-hcost(inssim,hnode,2-int(famsim)),cminp+mcdcash(hnode,2-int(famsim)))
				COHpsim=COHpsim-premium
				if (COHpsim>cminp+mcdcash(hnode,2-int(famsim))-premium) then
					socinspsim=0
				else
					socinspsim=1
				endif
			elseif (inssim.eq.3) then
				COHpsim = max(assetpsim+incomep(wpnode,simage)-hcost(inssim,hnode,2-int(famsim)),cminp+mcdcash(hnode,2-int(famsim)))
				COHpsim=COHpsim-premiumPO
				if (COHpsim>cminp+mcdcash(hnode,2-int(famsim))-premiumPO) then
					socinspsim=0
				else
					socinspsim=1
				endif
			else
				COHpsim = max(assetpsim+incomep(wpnode,simage)-hcost(inssim,hnode,2-int(famsim)),cminp+mcdcash(hnode,2-int(famsim)))
				if (COHpsim>cminp+mcdcash(hnode,2-int(famsim))) then
					socinspsim=0
				else
					socinspsim=1
				endif					
			endif


			if (LFPsim.eq.3) then
				COHksim = max(assetksim+incomek(wknode,simage)           ,cmink)
			elseif (LFPsim.eq.2) then
				COHksim = max(assetksim+PTwagefrac*incomek(wknode,simage),cmink)
			else
				COHksim = max(assetksim,cmink)						
			endif					
			if (COHksim==cmink) then
				socinsksim=1
			else
				socinsksim=0
			endif


			if (marsim==1) then 
				savsim =COHpsim+COHksim-conexppsim-conexpksim
				savpsim=kapsim*savsim !if last period, kapsim=0 so will look like no bequest...
				savksim=(1.0-kapsim)*savsim
			else
				savpsim=COHpsim-conexppsim
				savksim=COHksim-conexpksim
			endif

			if (savpsim<0) then
				conexppsim=conexppsim+savpsim
				savpsim=0
			endif
			if (savksim<0) then
				conexpksim=conexpksim+savksim
				savksim=0
			endif

			if (diethispd.eq.1) then
				beqsim=r*savpsim   
				savksim = savksim+savpsim
			else
				beqsim=0.0
				ktermasim=0.0
			endif

			!record matrix of all sorts of values
			inctypep(index)  = wpnode
			inctypek(index)  = incktype
			xval_ins(index,t) = xvaluemat(1)
			xval_noins(index,t) = xvaluemat(2)
			if (ni==3) then
				xval_insPO(index,t) = xvaluemat(3)
			endif
			xvalp(index,t) 	 = xvaluep
			xvalk(index,t)   = xvaluek
			conk(index,t)    = conexpksim
			conp(index,t)    = conexppsim
			savp(index,t)    = savpsim
			savk(index,t)	 	 = savksim
			assp(index,t)		 = assetpsim
			assk(index,t)		 = assetksim
			COHp(index,t)		 = COHpsim
			COHk(index,t)		 = COHksim
			ssip(index,t)	   = socinspsim
			ssik(index,t)    = socinsksim
			incp(index,t)    = incomep(wpnode,simage)
			inck(index,t)    = incomek(wknode,simage)
			if (inssim==1) then 
				insd(index,t)      = 1	
			elseif (inssim==3) then
				insdPO(index,t) = 1
			endif
			
			labor(index,t)   = LFPsim
			hlth(index,t)    = hnode
			informal(index,t)= famsim
			marstat(index,t) = marsim
			muwgt(index,t)   = musim
			kappct(index,t)	 = kapsim
			bequest(index,t) = beqsim
			kterma(index,t)  = ktermasim
			wstardiv(index,t)  = wstar_divorced
			wstarmar(index,t)  = wstar_married

			if (diethispd.eq.1) dead=1

		endif	
		
	500 continue
	

		
	enddo	!end loop time (t)
	! _____________________________________


contains 
! *******************************************************


! *******************************************************
subroutine find_valpwgt_sim(val, &
                            xvalue_adj,vp_adj,vk_adj,cp_adj,ck_adj,kap_adj,mu_adj,LFP_adj,fam_adj)

	integer m, i
	integer mu_adj
	real(prec) val, vp_adj,vk_adj,cp_adj,ck_adj,kap_adj,LFP_adj,fam_adj

	real(prec) VFadj(6,na,nonodesk), VFpadj(6,na,nonodesk), VFkadj(6,na,nonodesk)
	real(prec) Cpadj(6,na,nonodesk), Ckadj(6,na,nonodesk), kapadj(6,na,nonodesk)
	real(prec) xCVF_adj(6),xvalue_adj,xvaluep_adj
	integer locadj(1)

 do m=musimstart,nm

 	!solve for value and then valuep at that m
	VFadj(1,:,:)=valueNT0_married(simage,ins,bnode,m,:,:)
	VFadj(2,:,:)=valueNT100_married(simage,ins,bnode,m,:,:)
	VFadj(3,:,:)=valuePT0_married(simage,ins,bnode,m,:,:)
	VFadj(4,:,:)=valuePT100_married(simage,ins,bnode,m,:,:)
	VFadj(5,:,:)=valueFT0_married(simage,ins,bnode,m,:,:)
	VFadj(6,:,:)=valueFT100_married(simage,ins,bnode,m,:,:)

	VFpadj(1,:,:)=valuepNT0_married(simage,ins,bnode,m,:,:)
	VFpadj(2,:,:)=valuepNT100_married(simage,ins,bnode,m,:,:)
	VFpadj(3,:,:)=valuepPT0_married(simage,ins,bnode,m,:,:)
	VFpadj(4,:,:)=valuepPT100_married(simage,ins,bnode,m,:,:)
	VFpadj(5,:,:)=valuepFT0_married(simage,ins,bnode,m,:,:)
	VFpadj(6,:,:)=valuepFT100_married(simage,ins,bnode,m,:,:)

	VFkadj(1,:,:)=valuekNT0_married(simage,ins,bnode,m,:,:)
	VFkadj(2,:,:)=valuekNT100_married(simage,ins,bnode,m,:,:)
	VFkadj(3,:,:)=valuekPT0_married(simage,ins,bnode,m,:,:)
	VFkadj(4,:,:)=valuekPT100_married(simage,ins,bnode,m,:,:)
	VFkadj(5,:,:)=valuekFT0_married(simage,ins,bnode,m,:,:)
	VFkadj(6,:,:)=valuekFT100_married(simage,ins,bnode,m,:,:)

	VFadj(1,:,:)=valueNT0_married(simage,ins,bnode,m,:,:)
	VFadj(2,:,:)=valueNT100_married(simage,ins,bnode,m,:,:)
	VFadj(3,:,:)=valuePT0_married(simage,ins,bnode,m,:,:)
	VFadj(4,:,:)=valuePT100_married(simage,ins,bnode,m,:,:)
	VFadj(5,:,:)=valueFT0_married(simage,ins,bnode,m,:,:)
	VFadj(6,:,:)=valueFT100_married(simage,ins,bnode,m,:,:)

	Cpadj(1,:,:)=conexppNT0_married(simage,ins,bnode,m,:,:)
	Cpadj(2,:,:)=conexppNT100_married(simage,ins,bnode,m,:,:)
	Cpadj(3,:,:)=conexppPT0_married(simage,ins,bnode,m,:,:)
	Cpadj(4,:,:)=conexppPT100_married(simage,ins,bnode,m,:,:)
	Cpadj(5,:,:)=conexppFT0_married(simage,ins,bnode,m,:,:)
	Cpadj(6,:,:)=conexppFT100_married(simage,ins,bnode,m,:,:)
	
	Ckadj(1,:,:)=conexpkNT0_married(simage,ins,bnode,m,:,:)
	Ckadj(2,:,:)=conexpkNT100_married(simage,ins,bnode,m,:,:)
	Ckadj(3,:,:)=conexpkPT0_married(simage,ins,bnode,m,:,:)
	Ckadj(4,:,:)=conexpkPT100_married(simage,ins,bnode,m,:,:)
	Ckadj(5,:,:)=conexpkFT0_married(simage,ins,bnode,m,:,:)
	Ckadj(6,:,:)=conexpkFT100_married(simage,ins,bnode,m,:,:)

	kapadj(1,:,:)=kappaNT0_married(simage,ins,bnode,m,:,:)
	kapadj(2,:,:)=kappaNT100_married(simage,ins,bnode,m,:,:)
	kapadj(3,:,:)=kappaPT0_married(simage,ins,bnode,m,:,:)
	kapadj(4,:,:)=kappaPT100_married(simage,ins,bnode,m,:,:)
	kapadj(5,:,:)=kappaFT0_married(simage,ins,bnode,m,:,:)
	kapadj(6,:,:)=kappaFT100_married(simage,ins,bnode,m,:,:)

	xCVF_adj(1)=bilinear(xpts_bi,lb_bi,ub_bi,VFadj(1,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
	xCVF_adj(2)=bilinear(xpts_bi,lb_bi,ub_bi,VFadj(2,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
	xCVF_adj(3)=bilinear(xpts_bi,lb_bi,ub_bi,VFadj(3,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
	xCVF_adj(4)=bilinear(xpts_bi,lb_bi,ub_bi,VFadj(4,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
	xCVF_adj(5)=bilinear(xpts_bi,lb_bi,ub_bi,VFadj(5,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
	xCVF_adj(6)=bilinear(xpts_bi,lb_bi,ub_bi,VFadj(6,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))

	xvalue_adj=maxval(xCVF_adj)
	locadj=maxloc(xCVF_adj)
	xvaluep_adj=bilinear(xpts_bi,lb_bi,ub_bi,VFpadj(locadj(1),int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))

 	!check where valuep is good enough to satistfy parent
    if (xvaluep_adj>=val .or. m==nm) then

      mu_adj=m
      if (xvalue_adj==xCVF_adj(1)) then
        kap_adj=bilinear(xpts_bi,lb_bi,ub_bi,kapadj(1,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        vp_adj=xvaluep_adj
        vk_adj=bilinear(xpts_bi,lb_bi,ub_bi,VFkadj(1,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        cp_adj=bilinear(xpts_bi,lb_bi,ub_bi,Cpadj(1,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        ck_adj=bilinear(xpts_bi,lb_bi,ub_bi,Ckadj(1,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        LFP_adj=1
        fam_adj=0        
      elseif (xvalue_adj==xCVF_adj(2)) then
        kap_adj=bilinear(xpts_bi,lb_bi,ub_bi,kapadj(2,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        vp_adj=xvaluep_adj
        vk_adj=bilinear(xpts_bi,lb_bi,ub_bi,VFkadj(2,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        cp_adj=bilinear(xpts_bi,lb_bi,ub_bi,Cpadj(2,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        ck_adj=bilinear(xpts_bi,lb_bi,ub_bi,Ckadj(2,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        LFP_adj=1
        fam_adj=1
      elseif (xvalue_adj==xCVF_adj(3)) then
        kap_adj=bilinear(xpts_bi,lb_bi,ub_bi,kapadj(3,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        vp_adj=xvaluep_adj
        vk_adj=bilinear(xpts_bi,lb_bi,ub_bi,VFkadj(3,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        cp_adj=bilinear(xpts_bi,lb_bi,ub_bi,Cpadj(3,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        ck_adj=bilinear(xpts_bi,lb_bi,ub_bi,Ckadj(3,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        LFP_adj=2
        fam_adj=0       
      elseif (xvalue_adj==xCVF_adj(4)) then
        kap_adj=bilinear(xpts_bi,lb_bi,ub_bi,kapadj(4,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        vp_adj=xvaluep_adj
        vk_adj=bilinear(xpts_bi,lb_bi,ub_bi,VFkadj(4,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        cp_adj=bilinear(xpts_bi,lb_bi,ub_bi,Cpadj(4,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        ck_adj=bilinear(xpts_bi,lb_bi,ub_bi,Ckadj(4,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        LFP_adj=2
        fam_adj=1
      elseif (xvalue_adj==xCVF_adj(5)) then
        kap_adj=bilinear(xpts_bi,lb_bi,ub_bi,kapadj(5,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        vp_adj=xvaluep_adj
        vk_adj=bilinear(xpts_bi,lb_bi,ub_bi,VFkadj(5,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        cp_adj=bilinear(xpts_bi,lb_bi,ub_bi,Cpadj(5,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        ck_adj=bilinear(xpts_bi,lb_bi,ub_bi,Ckadj(5,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        LFP_adj=3
        fam_adj=0       
      elseif (xvalue_adj==xCVF_adj(6)) then
        kap_adj=bilinear(xpts_bi,lb_bi,ub_bi,kapadj(6,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        vp_adj=xvaluep_adj
        vk_adj=bilinear(xpts_bi,lb_bi,ub_bi,VFkadj(6,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        cp_adj=bilinear(xpts_bi,lb_bi,ub_bi,Cpadj(6,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        ck_adj=bilinear(xpts_bi,lb_bi,ub_bi,Ckadj(6,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        LFP_adj=3
        fam_adj=1
      endif

      goto 1234
    endif

  enddo

1234 continue
end subroutine find_valpwgt_sim
! *******************************************************


! *******************************************************
subroutine find_valkwgt_sim(val, &
                            xvalue_adj,vp_adj,vk_adj,cp_adj,ck_adj,kap_adj,mu_adj,LFP_adj,fam_adj)

	integer m, i
	integer mu_adj
	real(prec) val, v_adj,vp_adj,vk_adj,cp_adj,ck_adj,kap_adj,LFP_adj,fam_adj

	real(prec) VFadj(6,na,nonodesk), VFpadj(6,na,nonodesk), VFkadj(6,na,nonodesk)
	real(prec) Cpadj(6,na,nonodesk), Ckadj(6,na,nonodesk), kapadj(6,na,nonodesk)
	real(prec) xCVF_adj(6),xvalue_adj,xvaluek_adj
	integer locadj(1)

 do m=musimstart,1,-1

 	!solve for value and then valuep at that m
	VFadj(1,:,:)=valueNT0_married(simage,ins,bnode,m,:,:)
	VFadj(2,:,:)=valueNT100_married(simage,ins,bnode,m,:,:)
	VFadj(3,:,:)=valuePT0_married(simage,ins,bnode,m,:,:)
	VFadj(4,:,:)=valuePT100_married(simage,ins,bnode,m,:,:)
	VFadj(5,:,:)=valueFT0_married(simage,ins,bnode,m,:,:)
	VFadj(6,:,:)=valueFT100_married(simage,ins,bnode,m,:,:)

	VFpadj(1,:,:)=valuepNT0_married(simage,ins,bnode,m,:,:)
	VFpadj(2,:,:)=valuepNT100_married(simage,ins,bnode,m,:,:)
	VFpadj(3,:,:)=valuepPT0_married(simage,ins,bnode,m,:,:)
	VFpadj(4,:,:)=valuepPT100_married(simage,ins,bnode,m,:,:)
	VFpadj(5,:,:)=valuepFT0_married(simage,ins,bnode,m,:,:)
	VFpadj(6,:,:)=valuepFT100_married(simage,ins,bnode,m,:,:)

	VFkadj(1,:,:)=valuekNT0_married(simage,ins,bnode,m,:,:)
	VFkadj(2,:,:)=valuekNT100_married(simage,ins,bnode,m,:,:)
	VFkadj(3,:,:)=valuekPT0_married(simage,ins,bnode,m,:,:)
	VFkadj(4,:,:)=valuekPT100_married(simage,ins,bnode,m,:,:)
	VFkadj(5,:,:)=valuekFT0_married(simage,ins,bnode,m,:,:)
	VFkadj(6,:,:)=valuekFT100_married(simage,ins,bnode,m,:,:)

	VFadj(1,:,:)=valueNT0_married(simage,ins,bnode,m,:,:)
	VFadj(2,:,:)=valueNT100_married(simage,ins,bnode,m,:,:)
	VFadj(3,:,:)=valuePT0_married(simage,ins,bnode,m,:,:)
	VFadj(4,:,:)=valuePT100_married(simage,ins,bnode,m,:,:)
	VFadj(5,:,:)=valueFT0_married(simage,ins,bnode,m,:,:)
	VFadj(6,:,:)=valueFT100_married(simage,ins,bnode,m,:,:)

	Cpadj(1,:,:)=conexppNT0_married(simage,ins,bnode,m,:,:)
	Cpadj(2,:,:)=conexppNT100_married(simage,ins,bnode,m,:,:)
	Cpadj(3,:,:)=conexppPT0_married(simage,ins,bnode,m,:,:)
	Cpadj(4,:,:)=conexppPT100_married(simage,ins,bnode,m,:,:)
	Cpadj(5,:,:)=conexppFT0_married(simage,ins,bnode,m,:,:)
	Cpadj(6,:,:)=conexppFT100_married(simage,ins,bnode,m,:,:)
	
	Ckadj(1,:,:)=conexpkNT0_married(simage,ins,bnode,m,:,:)
	Ckadj(2,:,:)=conexpkNT100_married(simage,ins,bnode,m,:,:)
	Ckadj(3,:,:)=conexpkPT0_married(simage,ins,bnode,m,:,:)
	Ckadj(4,:,:)=conexpkPT100_married(simage,ins,bnode,m,:,:)
	Ckadj(5,:,:)=conexpkFT0_married(simage,ins,bnode,m,:,:)
	Ckadj(6,:,:)=conexpkFT100_married(simage,ins,bnode,m,:,:)

	kapadj(1,:,:)=kappaNT0_married(simage,ins,bnode,m,:,:)
	kapadj(2,:,:)=kappaNT100_married(simage,ins,bnode,m,:,:)
	kapadj(3,:,:)=kappaPT0_married(simage,ins,bnode,m,:,:)
	kapadj(4,:,:)=kappaPT100_married(simage,ins,bnode,m,:,:)
	kapadj(5,:,:)=kappaFT0_married(simage,ins,bnode,m,:,:)
	kapadj(6,:,:)=kappaFT100_married(simage,ins,bnode,m,:,:)

	xCVF_adj(1)=bilinear(xpts_bi,lb_bi,ub_bi,VFadj(1,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
	xCVF_adj(2)=bilinear(xpts_bi,lb_bi,ub_bi,VFadj(2,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
	xCVF_adj(3)=bilinear(xpts_bi,lb_bi,ub_bi,VFadj(3,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
	xCVF_adj(4)=bilinear(xpts_bi,lb_bi,ub_bi,VFadj(4,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
	xCVF_adj(5)=bilinear(xpts_bi,lb_bi,ub_bi,VFadj(5,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
	xCVF_adj(6)=bilinear(xpts_bi,lb_bi,ub_bi,VFadj(6,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))

	xvalue_adj=maxval(xCVF_adj)
	locadj=maxloc(xCVF_adj)
	xvaluek_adj=bilinear(xpts_bi,lb_bi,ub_bi,VFkadj(locadj(1),int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))

 	!check where valuep is good enough to satistfy parent
    if (xvaluek_adj>=val .or. m==1) then

      mu_adj=m
      if (xvalue_adj==xCVF_adj(1)) then
        kap_adj=bilinear(xpts_bi,lb_bi,ub_bi,kapadj(1,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        vp_adj=bilinear(xpts_bi,lb_bi,ub_bi,VFpadj(1,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        vk_adj=xvaluek_adj
        cp_adj=bilinear(xpts_bi,lb_bi,ub_bi,Cpadj(1,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        ck_adj=bilinear(xpts_bi,lb_bi,ub_bi,Ckadj(1,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        LFP_adj=1
        fam_adj=0        
      elseif (xvalue_adj==xCVF_adj(2)) then
        kap_adj=bilinear(xpts_bi,lb_bi,ub_bi,kapadj(2,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        vp_adj=bilinear(xpts_bi,lb_bi,ub_bi,VFpadj(2,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        vk_adj=xvaluek_adj
        cp_adj=bilinear(xpts_bi,lb_bi,ub_bi,Cpadj(2,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        ck_adj=bilinear(xpts_bi,lb_bi,ub_bi,Ckadj(2,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        LFP_adj=1
        fam_adj=1
      elseif (xvalue_adj==xCVF_adj(3)) then
        kap_adj=bilinear(xpts_bi,lb_bi,ub_bi,kapadj(3,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        vp_adj=bilinear(xpts_bi,lb_bi,ub_bi,VFpadj(3,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        vk_adj=xvaluek_adj
        cp_adj=bilinear(xpts_bi,lb_bi,ub_bi,Cpadj(3,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        ck_adj=bilinear(xpts_bi,lb_bi,ub_bi,Ckadj(3,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        LFP_adj=2
        fam_adj=0       
      elseif (xvalue_adj==xCVF_adj(4)) then
        kap_adj=bilinear(xpts_bi,lb_bi,ub_bi,kapadj(4,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        vp_adj=bilinear(xpts_bi,lb_bi,ub_bi,VFpadj(4,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        vk_adj=xvaluek_adj
        cp_adj=bilinear(xpts_bi,lb_bi,ub_bi,Cpadj(4,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        ck_adj=bilinear(xpts_bi,lb_bi,ub_bi,Ckadj(4,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        LFP_adj=2
        fam_adj=1
      elseif (xvalue_adj==xCVF_adj(5)) then
        kap_adj=bilinear(xpts_bi,lb_bi,ub_bi,kapadj(5,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        vp_adj=bilinear(xpts_bi,lb_bi,ub_bi,VFpadj(5,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        vk_adj=xvaluek_adj
        cp_adj=bilinear(xpts_bi,lb_bi,ub_bi,Cpadj(5,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        ck_adj=bilinear(xpts_bi,lb_bi,ub_bi,Ckadj(5,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        LFP_adj=3
        fam_adj=0       
      elseif (xvalue_adj==xCVF_adj(6)) then
        kap_adj=bilinear(xpts_bi,lb_bi,ub_bi,kapadj(6,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        vp_adj=bilinear(xpts_bi,lb_bi,ub_bi,VFpadj(6,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        vk_adj=xvaluek_adj
        cp_adj=bilinear(xpts_bi,lb_bi,ub_bi,Cpadj(6,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        ck_adj=bilinear(xpts_bi,lb_bi,ub_bi,Ckadj(6,int(lb_bi(1)):int(ub_bi(1)),int(lb_bi(2)):int(ub_bi(2))))
        LFP_adj=3
        fam_adj=1
      endif
      goto 1235
    endif

  enddo

1235 continue
end subroutine find_valkwgt_sim
! *******************************************************


!*******************************************************
subroutine find_wstar_div(fval,fmina,wgt1,wgt2)

    integer fage, fmina, jl, jm, ju
    real(prec) fval, a1,a2, wgt1,wgt2
    real(prec) y1,y2,iwgt1,iwgt2,guess

    iwgt1 = ((x_ak-ub_ak))/((lb_ak-ub_ak))
    iwgt2 = ((x_ak-lb_ak))/((ub_ak-lb_ak))

    jl=0
    ju=na+1
812  if (ju-jl.gt.1) then
      jm = (ju+jl)/2

      y1 = valuep_divorced(simage,2,bnode,jm,lb_ak)
      y2 = valuep_divorced(simage,2,bnode,jm,ub_ak)
      guess=iwgt1*y1+iwgt2*y2

	      if (fval.gt.guess) then
	          jl = jm
	      else
	          ju = jm
	      endif
	      goto 812
    	endif
    fmina = jl

    !find how far the asset amount lies between two nodes
    if (fmina.eq.na) then         
      fmina = na-1
    elseif (fmina.eq.0) then
      fmina = 1
    endif

    y1 = valuep_divorced(simage,2,bnode,fmina,lb_ak)
    y2 = valuep_divorced(simage,2,bnode,fmina,ub_ak)
    a1=iwgt1*y1+iwgt2*y2
         
    y1 = valuep_divorced(simage,2,bnode,fmina+1,lb_ak)
    y2 = valuep_divorced(simage,2,bnode,fmina+1,ub_ak)
    a2=iwgt1*y1+iwgt2*y2

    wgt1=(fval-a2)/(a1-a2)        !get weight on lower node
    wgt2=(fval-a1)/(a2-a1)        !get weight on higher node

end subroutine find_wstar_div
! *******************************************************

!*******************************************************

subroutine find_wstar_mar(pork,fval,fmina,wgt1,wgt2) !THIS ASSUMES MONEY ACCRUES TO PARENT

    integer pork,fage, fmina, jl, jm, ju
    real(prec) fval, a1,a2, wgt1,wgt2
    real(prec) y1,y2,iwgt1,iwgt2,guess

    jl=0
    ju=na+1
812 if (ju-jl.gt.1) then
      	jm = (ju+jl)/2

      	if (pork==1) then
	    	y1 = value(simage,2,bnode,musimstart,jm,lb_ak) 
	      	y2 = value(simage,2,bnode,musimstart,jm,ub_ak)
	      	iwgt1 = ((x_ak-ub_ak))/((lb_ak-ub_ak))
	      	iwgt2 = ((x_ak-lb_ak))/((ub_ak-lb_ak))
      	else 
	      	y1 = value(simage,2,bnode,musimstart,lb_ap,nodegridk+jm)
	      	y2 = value(simage,2,bnode,musimstart,ub_ap,nodegridk+jm)
	      	iwgt1 = ((x_ap-ub_ap))/((lb_ap-ub_ap))
	      	iwgt2 = ((x_ap-lb_ap))/((ub_ap-lb_ap))
      	endif
      	guess=iwgt1*y1+iwgt2*y2

      if (fval.gt.guess) then
          jl = jm
      else
          ju = jm
      endif
      goto 812
    endif
    fmina = jl

    !find how far the asset amount lies between two nodes
    if (fmina.eq.na) then         
      fmina = na-1
    elseif (fmina.eq.0) then
      fmina = 1
    endif

  	if (pork==1) then
	    y1 = value(simage,2,bnode,musimstart,fmina,lb_ak)
	    y2 = value(simage,2,bnode,musimstart,fmina,ub_ak)
	    a1=iwgt1*y1+iwgt2*y2

	    y1 = value(simage,2,bnode,musimstart,fmina+1,lb_ak)
	    y2 = value(simage,2,bnode,musimstart,fmina+1,ub_ak)
	    a2=iwgt1*y1+iwgt2*y2

	    wgt1=(fval-a2)/(a1-a2)        !get weight on lower node
	    wgt2=(fval-a1)/(a2-a1)        !get weight on higher node
    else
	    y1 = value(simage,2,bnode,musimstart,lb_ap,nodegridk+fmina)
	    y2 = value(simage,2,bnode,musimstart,ub_ap,nodegridk+fmina)
	    a1=iwgt1*y1+iwgt2*y2

	    y1 = value(simage,2,bnode,musimstart,lb_ap,nodegridk+fmina+1)
	    y2 = value(simage,2,bnode,musimstart,ub_ap,nodegridk+fmina+1)
	    a2=iwgt1*y1+iwgt2*y2

	    wgt1=(fval-a2)/(a1-a2)        !get weight on lower node
	    wgt2=(fval-a1)/(a2-a1)        !get weight on higher node
  	endif


end subroutine find_wstar_mar
! *******************************************************

end subroutine run_simperson


end module simulate
